#Amelia Gaudio and Ryan Welch
#12 August 2024
#Abortion Project
#Replication File

###Preamble
######################################################################################
rm(list=ls())

#check if packages are installed, and if not, install them
packages <- c("foreign", "Hmisc", "dplyr", "ggplot2", "GGally", "MatchIt", "sjPlot", "sjmisc", "margins", "interplot", "maps", "stargazer", "psych", "plm", "car", "cem", "rbounds", "texreg", "data.table")
if (length(setdiff(packages, rownames(installed.packages()))) > 0) {
  install.packages(setdiff(packages, rownames(installed.packages())), repos="http://cran.rstudio.com/")  
}

#attach packages
lapply(packages, library, character.only = T)

#default the stringsAsFactors argument to FALSE for all functions with that argument
options(stringsAsFactors=FALSE)

#turn off scientific notation
options(scipen=999)

#set seed
RNGkind("L'Ecuyer-CMRG")
set.seed(33603)

#office
#setwd("C:/Users/rwelch/Dropbox/Research/Abortion")

#Personal
setwd("C:/Users/ryanm/Dropbox/Research/Abortion")

#Load data
data_abo <- read.csv("abo_data_06Aug2024.csv", stringsAsFactors = F, strip.white = T)
################################################################################

# Save the plot as a PNG file with 300 dpi resolution
png("Figure1_HiRes.png", width = 6, height = 4, units = "in", res = 300)


#Figure 1. Distribution of Dependent Variable
props <- prop.table(table(data_abo$index3))
barplot(props,
        xlab="Abortion Rights Index Score",
        ylab = "Proportion",
        ylim = (c(0, 0.3))
)

dev.off()


png("Figure2_HiRes.png", width = 6, height = 4, units = "in", res = 300)

#Figure 2. Proportion of Country-Years with CEDAW-OP Ratification
props <- prop.table(table(data_abo$CEDAWOpt))
barplot(props,
        xlab="Ratification Status",
        ylab = "Proportion"
)

dev.off()


png("Figure3_HiRes.png", width = 6, height = 4, units = "in", res = 300)

#Figure 3. Distribution of the Women's Civil Society Participation Index
#Let's do a civil society variable together.
hist(data_abo$v2x_gencs_VDEM,
     freq=F,
     xlab="Women's Civil Society Participation Score",
     ylab="Probability Density",
     main="") 

dev.off()



#Create matched dataset on the Treaty of Interest (CEDAW-OP)

#Missing values not allowed in the treatment variable
data_abo2 <- data_abo[,c("cowcode", "year", "index3", "CEDAWOpt", 
                         "v2x_gencs_VDEM", "polity2.x", "hasNHRI", "injud_CR",
                         "physint_CR", "women_parl", "chrpc_RCS", "muspc_RCS",
                         "lgdp", "lpop")]

data_abo3 <- data_abo2[complete.cases(data_abo2),]
##################################################################################

#Nearest Neighbor (one to one)
nearest.out1<-matchit(CEDAWOpt ~ v2x_gencs_VDEM+ polity2.x + hasNHRI + 
                        injud_CR + physint_CR + women_parl + chrpc_RCS + 
                        muspc_RCS + lgdp + lpop,
                      data = data_abo3, method ="nearest", distance = "logit", 
                      discard = "both", reestimate = TRUE)

nearest.out1


#Table A3
summary(nearest.out1)
#distance              0.5518        0.3962

#create data frame
nearest.data1<-match.data(nearest.out1)

#Regression
nearestMatchOP_RegM<-lm(index3 ~ CEDAWOpt*v2x_gencs_VDEM + polity2.x + hasNHRI + 
                         injud_CR + physint_CR + women_parl + chrpc_RCS + 
                         muspc_RCS + lgdp + lpop, data = nearest.data1)

#Table 1 (NN(1:1))
summary(nearestMatchOP_RegM)

#sensitivity analysis
treat <- nearest.data1$CEDAWOpt
outcome <- nearest.data1$index3
psens_resultNN <- psens(outcome, treat)
print(psens_resultNN)


#Table A8 (Nearest Neighbor (1:1))
psens(outcome, treat, Gamma = 3, GammaInc = 0.1)
# 2.2           0      0.0129
# 2.3           0      0.0652

vif(nearestMatchOP_RegM)

#Nothing over 5, much less 10. Highest is population (3.7) (except of course for
#constituents, which we would expect)


##############################################################################

###Optimal
optimal.out1<-matchit(CEDAWOpt ~ v2x_gencs_VDEM+ polity2.x + hasNHRI +
                        injud_CR + physint_CR + women_parl + chrpc_RCS + 
                        muspc_RCS + lgdp + lpop,
                      data = data_abo3, method ="optimal",  
                      discard = "both")

optimal.out1

#Table A4
summary(optimal.out1)
#distance              0.5540        0.3951

#create data frame
optimal.data1<-match.data(optimal.out1)

#sensitivity analysis
treat <- optimal.data1$CEDAWOpt
outcome <- optimal.data1$index3
psens_resultOptimal <- psens(outcome, treat)
print(psens_resultOptimal)

#Table A8 (Optimal)
psens(outcome, treat, Gamma = 3, GammaInc = 0.1)
# 2.2           0      0.0164
# 2.3           0      0.0782

#Regression
optimalMatchOP_RegM<-lm(index3 ~ CEDAWOpt*v2x_gencs_VDEM + polity2.x + hasNHRI + 
                         injud_CR + physint_CR + women_parl + chrpc_RCS + 
                         muspc_RCS + lgdp + lpop, data = optimal.data1)


#Table 1 (Optimal)
summary(optimalMatchOP_RegM)





vif(optimalMatchOP_RegM)

#Effectively the same as NN

################################################################################

#Full

full.out1<-matchit(CEDAWOpt ~ v2x_gencs_VDEM+ polity2.x + hasNHRI + 
                     injud_CR + physint_CR + women_parl + chrpc_RCS + 
                     muspc_RCS + lgdp + lpop,
                   data = data_abo3, method ="full",  
                   discard = "both")

full.out1

#Table A5
summary(full.out1)
#distance              0.5540        0.5540 

#create data frame
full.data1<-match.data(full.out1)

#sensitivity analysis
treat <- full.data1$CEDAWOpt
outcome <- full.data1$index3
psens_resultFull <- psens(outcome, treat)
print(psens_resultFull)

#Table A8 (Optimal Full)
psens(outcome, treat, Gamma = 3, GammaInc = 0.1)
# 2.3           0      0.0317
# 2.4           0      0.1348

#Regression
fullMatchOP_RegM<-lm(index3 ~ CEDAWOpt*v2x_gencs_VDEM + polity2.x + hasNHRI + 
                      injud_CR + physint_CR + women_parl + chrpc_RCS + 
                      muspc_RCS + lgdp + lpop, data = full.data1)


#Table 1 (Full)
summary(fullMatchOP_RegM)

vif(fullMatchOP_RegM)


#Figure 4.

png("Figure4_HiRes.png", width = 6, height = 4, units = "in", res = 300)


#Effect of CEDAW-OP Ratification conditional on Civil Society Strength
interplot(m = fullMatchOP_RegM, var1 = "CEDAWOpt", var2 = "v2x_gencs_VDEM", ci = 0.90) +
  #Add labels for X and Y axes
  xlab("Women's Civil Society Participation") +
  ylab("Estimated Effect of CEDAW-OP") +
  #Add a horizontal line at y = 0
  geom_hline(yintercept = 0, linetype = "dashed")

dev.off()

mean(data_abo$v2x_cspart_VDEM>=0.70, na.rm=T)
#[1] 0.3884094

#############################################################################

#Genetic

genetic.out1<-matchit(CEDAWOpt ~ v2x_gencs_VDEM+ polity2.x + hasNHRI +
                        injud_CR + physint_CR + women_parl + chrpc_RCS + 
                        muspc_RCS + lgdp + lpop,
                      data = data_abo3, method ="genetic", 
                      discard = "both", pop.size=1000)

genetic.out1

#Table A6
summary(genetic.out1)
#distance              0.5540        0.3763

#create data frame
genetic.data1<-match.data(genetic.out1)

#sensitivity analysis
treat <- genetic.data1$CEDAWOpt
outcome <- genetic.data1$index3
psens_resultGen <- psens(outcome, treat)
print(psens_resultGen)


#Table A8 (Genetic)
psens(outcome, treat, Gamma = 3, GammaInc = 0.1)
# 2.4           0      0.0300
# 2.5           0      0.1084


#Regression
geneticMatchOP_RegM<-lm(index3 ~ CEDAWOpt*v2x_gencs_VDEM + polity2.x + hasNHRI + 
                         injud_CR + physint_CR + women_parl + chrpc_RCS + 
                         muspc_RCS + lgdp + lpop, data = genetic.data1)

#Table 1 (Genetic)
summary(geneticMatchOP_RegM)
#Note: Estimates may differ from run to run. Even when one sets the seed, the
#genetic algorithm is not deterministic (see Sekhon and Mebane 1998; Diamond and
#Sekhon 2013; and Goldberg 1991 for more details). That said, every time we
#estimated it, our results told the same substantive story (i.e., positive beta 
#coefficients) at statistically significant levels.


vif(geneticMatchOP_RegM)



################################################################################


#Cardinality

cardinality.out1<-matchit(CEDAWOpt ~ v2x_gencs_VDEM+ polity2.x + hasNHRI + 
                            injud_CR + physint_CR + women_parl + chrpc_RCS + 
                            muspc_RCS + lgdp + lpop,
                          data = data_abo3, method = "cardinality", tols = 0.43)

cardinality.out1

#Table A7
summary(cardinality.out1)


#create data frame
cardinality.data1<-match.data(cardinality.out1)

#sensitivity analysis
treat <- cardinality.data1$CEDAWOpt
outcome <- cardinality.data1$index3
psens_resultCard <- psens(outcome, treat)
print(psens_resultCard)

#Table A8 (Cardinality)
psens(outcome, treat, Gamma = 3, GammaInc = 0.1)
# 2.4           0      0.0287
# 2.5           0      0.1065

#Regression
cardinalityMatchOP_RegM<-lm(index3 ~ CEDAWOpt*v2x_gencs_VDEM + polity2.x + hasNHRI + 
                             injud_CR + physint_CR + women_parl + chrpc_RCS + 
                             muspc_RCS + lgdp + lpop, data = cardinality.data1)

#Table 1 (Cardinality)
summary(cardinalityMatchOP_RegM)


vif(cardinalityMatchOP_RegM)
#Same

###########################################################################



htmlreg(
  list(nearestMatchOP_RegM, optimalMatchOP_RegM, fullMatchOP_RegM,
       geneticMatchOP_RegM, cardinalityMatchOP_RegM)
)

stargazer(nearestMatchOP_RegM, optimalMatchOP_RegM, fullMatchOP_RegM, type="html")

stargazer(geneticMatchOP_RegM, cardinalityMatchOP_RegM,
          type = "html")

stargazer(nearestMatchOP_RegM, fullMatchOP_RegM, geneticMatchOP_RegM, type="html")


group1_models <- list(nearestMatchOP_RegM, optimalMatchOP_RegM, fullMatchOP_RegM)
group2_models <- list(geneticMatchOP_RegM, cardinalityMatchOP_RegM)
stargazer(group1_models, group2_models, type="html")


################################################################################

###All robustness checks with (Optimal) Full Matching due to its superior balance.

##############################################################################

#With CEDAW control

#Missing values not allowed in the treatment variable
data_abo6 <- data_abo[,c("cowcode", "year", "index3", "CEDAWOpt", "CEDAW", 
                         "v2x_gencs_VDEM", "polity2.x", "hasNHRI", "injud_CR",
                         "physint_CR", "women_parl", "chrpc_RCS", "muspc_RCS",
                         "lgdp", "lpop")]

data_abo7 <- data_abo6[complete.cases(data_abo6),]



################################################################################

#Full

full.outA2<-matchit(CEDAWOpt ~ v2x_gencs_VDEM+ polity2.x + hasNHRI + 
                      injud_CR + physint_CR + women_parl + chrpc_RCS + 
                      muspc_RCS + lgdp + lpop + CEDAW,
                    data = data_abo7, method ="full",  
                    discard = "both")

##Table 1##
full.outA2
summary(full.outA2)
#distance              0.5585        0.5584 

#create data frame
full.dataA2<-match.data(full.outA2)

#Regression
fullMatchOP_RegMA2<-lm(index3 ~ CEDAWOpt*v2x_gencs_VDEM + polity2.x + hasNHRI + 
                        injud_CR + physint_CR + women_parl + chrpc_RCS + 
                        muspc_RCS + lgdp + lpop + CEDAW, data = full.dataA2)

summary(fullMatchOP_RegMA2)
#CEDAWOpt:v2x_gencs_VDEM  5.043868   0.961800   5.244  0.00000017225403568 ***

stargazer(fullMatchOP_RegMA2, type="html")




#With Forman-Rabinovici and Sommer (2018) coding rules

#Create matched dataset on the Treaty of Interest (CEDAW-OP)

#Missing values not allowed in the treatment variable
data_abo4 <- data_abo[,c("cowcode", "year", "Index1", "CEDAWOpt", 
                         "v2x_gencs_VDEM", "polity2.x", "hasNHRI", "injud_CR",
                         "physint_CR", "women_parl", "chrpc_RCS", "muspc_RCS",
                         "lgdp", "lpop")]

data_abo5 <- data_abo4[complete.cases(data_abo4),]


################################################################################

#Full

full.outA<-matchit(CEDAWOpt ~ v2x_gencs_VDEM+ polity2.x + hasNHRI + 
                     injud_CR + physint_CR + women_parl + chrpc_RCS + 
                     muspc_RCS + lgdp + lpop,
                   data = data_abo5, method ="full",  
                   discard = "both")

##Table 1##
full.outA
summary(full.outA)


#create data frame
full.dataA<-match.data(full.outA)

#Regression
fullMatchOP_RegMA<-lm(Index1 ~ CEDAWOpt*v2x_gencs_VDEM + polity2.x + hasNHRI + 
                      injud_CR + physint_CR + women_parl + chrpc_RCS + 
                      muspc_RCS + lgdp + lpop, data = full.dataA)

summary(fullMatchOP_RegMA)
#CEDAWOpt:v2x_gencs_VDEM  1.812039   0.521167   3.477             0.000516 ***

stargazer(fullMatchOP_RegMA, type="html")


#############################################################################

#Without interaction term
fullMatchOP_RegMA3<-lm(index3 ~ CEDAWOpt + v2x_gencs_VDEM + polity2.x + hasNHRI + 
                       injud_CR + physint_CR + women_parl + chrpc_RCS + 
                       muspc_RCS + lgdp + lpop, data = full.data1)

summary(fullMatchOP_RegMA3)

stargazer(fullMatchOP_RegMA3, type="html") 


##############################################################################

#Tables for Memo including execrlc (Table M1)

#Create matched dataset on the Treaty of Interest (CEDAW-OP)

#Missing values not allowed in the treatment variable
data_abo8 <- data_abo[,c("cowcode", "year", "index3", "CEDAWOpt", 
                         "v2x_gencs_VDEM", "polity2.x", "hasNHRI", "injud_CR", "execrlc",
                         "physint_CR", "women_parl", "chrpc_RCS", "muspc_RCS",
                         "lgdp", "lpop")]

data_abo9 <- data_abo8[complete.cases(data_abo8),]
##################################################################################

#Nearest Neighbor (one to one)
nearest.outM<-matchit(CEDAWOpt ~ v2x_gencs_VDEM+ polity2.x + hasNHRI + 
                        injud_CR + execrlc + physint_CR + women_parl + chrpc_RCS + 
                        muspc_RCS + lgdp + lpop,
                      data = data_abo9, method ="nearest", distance = "logit", 
                      discard = "both", reestimate = TRUE)

##Table##
nearest.outM
summary(nearest.outM)
#distance              0.5518        0.3962

#create data frame
nearest.dataM<-match.data(nearest.outM)

#Regression
nearestMatchOP_RegM<-lm(index3 ~ CEDAWOpt*v2x_gencs_VDEM + polity2.x + hasNHRI + 
                         injud_CR + execrlc + physint_CR + women_parl + chrpc_RCS + 
                         muspc_RCS + lgdp + lpop, data = nearest.dataM)

summary(nearestMatchOP_RegM)
#CEDAWOpt:v2x_gencs_VDEM   2.1745935   1.0907107   1.994              0.04633 *  





##############################################################################

###Optimal
optimal.outM<-matchit(CEDAWOpt ~ v2x_gencs_VDEM+ polity2.x + hasNHRI +
                        injud_CR + execrlc + physint_CR + women_parl + chrpc_RCS + 
                        muspc_RCS + lgdp + lpop,
                      data = data_abo9, method ="optimal",  
                      discard = "both")

##Table 1##
optimal.outM
summary(optimal.outM)
#distance              0.5540        0.3951

#create data frame
optimal.dataM<-match.data(optimal.outM)

#Regression
optimalMatchOP_RegM<-lm(index3 ~ CEDAWOpt*v2x_gencs_VDEM + polity2.x + hasNHRI + 
                         injud_CR + execrlc + physint_CR + women_parl + chrpc_RCS + 
                         muspc_RCS + lgdp + lpop, data = optimal.dataM)

summary(optimalMatchOP_RegM)
#CEDAWOpt:v2x_gencs_VDEM   1.4665354   1.0897906   1.346             0.178568    


################################################################################

#Full

full.outM<-matchit(CEDAWOpt ~ v2x_gencs_VDEM+ polity2.x + hasNHRI + 
                     injud_CR + execrlc + physint_CR + women_parl + chrpc_RCS + 
                     muspc_RCS + lgdp + lpop,
                   data = data_abo9, method ="full",  
                   discard = "both")

##Table 1##
full.outM
summary(full.outM)
#distance              0.5540        0.5540 

#create data frame
full.dataM<-match.data(full.outM)

#Regression
fullMatchOP_RegM<-lm(index3 ~ CEDAWOpt*v2x_gencs_VDEM + polity2.x + hasNHRI + 
                      injud_CR + execrlc + physint_CR + women_parl + chrpc_RCS + 
                      muspc_RCS + lgdp + lpop, data = full.dataM)

summary(fullMatchOP_RegM)
#CEDAWOpt:v2x_gencs_VDEM  3.6523082  0.9445348   3.867             0.000113 ***


#############################################################################

#Genetic

genetic.outM<-matchit(CEDAWOpt ~ v2x_gencs_VDEM+ polity2.x + hasNHRI +
                        injud_CR + execrlc + physint_CR + women_parl + chrpc_RCS + 
                        muspc_RCS + lgdp + lpop,
                      data = data_abo9, method ="genetic", 
                      discard = "both", pop.size=1000)

##Table 1##
genetic.outM
summary(genetic.outM)
#distance              0.5540        0.3763

#create data frame
genetic.dataM<-match.data(genetic.outM)

#Regression
geneticMatchOP_RegM<-lm(index3 ~ CEDAWOpt*v2x_gencs_VDEM + polity2.x + hasNHRI + 
                         injud_CR + execrlc + physint_CR + women_parl + chrpc_RCS + 
                         muspc_RCS + lgdp + lpop, data = genetic.dataM)

summary(geneticMatchOP_RegM)
#CEDAWOpt:v2x_gencs_VDEM   3.5270114   1.0948691   3.221              0.00130 ** 



################################################################################


#Cardinality

cardinality.outM<-matchit(CEDAWOpt ~ v2x_gencs_VDEM+ polity2.x + hasNHRI + 
                            injud_CR + execrlc + physint_CR + women_parl + chrpc_RCS + 
                            muspc_RCS + lgdp + lpop,
                          data = data_abo9, method = "cardinality", tols = 0.43)

##Table 1##
cardinality.outM
summary(cardinality.outM)


#create data frame
cardinality.dataM<-match.data(cardinality.outM)

#Regression
cardinalityMatchOP_RegM<-lm(index3 ~ CEDAWOpt*v2x_gencs_VDEM + polity2.x + hasNHRI + 
                             injud_CR + execrlc + physint_CR + women_parl + chrpc_RCS + 
                             muspc_RCS + lgdp + lpop, data = cardinality.dataM)

summary(cardinalityMatchOP_RegM)
#CEDAWOpt:v2x_gencs_VDEM   3.7382622   1.0544744   3.545             0.000402 ***


###########################################################################



htmlreg(
  list(nearestMatchOP_RegM, optimalMatchOP_RegM, fullMatchOP_RegM,
       geneticMatchOP_RegM, cardinalityMatchOP_RegM)
)

